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ABSTRACT 

We investigate the influence of non equilibrium thermodynamics on 
cosmological structure formation. In this paper, we consider the collapse of 
planar perturbations usually called “Zel’dovich pancakes”. We have developed 
for that purpose a new two fluids (gas and dark matter) hydrodynamical 
code, with three different thermodynamical species: electrons, ions and neutral 
particles Tg 7 ^ T, 7 ^ T„). We describe in details the complex structure of 
accretion shock waves. We include several relevant processes for a low density, 
high temperature, collisional plasma such as non-equilibrium chemical reactions, 
cooling, shock heating, thermal energy equipartition between electrons, ions 
and neutral particles and electronic conduction. We hnd two different regions 
in the pancake structure: a thermal precursor ahead of the compression front 
and an equipartition wave after the compression front where electrons and 
ions temperatures differ signihcantly. This complex structure may have two 
interesting consequences: pre-heating of unshocked regions in the vicinity of 
massive X-ray clusters and ions and electrons temperatures differences in the 
outer regions of X-rays clusters. 


Subject headings: Cosmology: theory - hydrodynamics - methods: numerical 
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1. INTRODUCTION. 


The thermodynamical state of the Intergalactic Medium is of primary importance to 
study the formation and the evolution of cosmic structures. In the dense central regions 
of galaxy clusters, cooling is likely to play a dominant role. In the outer regions, the 
density is rather low and allows an adiabatic treatment of gas dynamics. In the same time, 
non-equilibrium thermodynamics occurs in this in this hot and diffuse plasma. Indeed, for 
high temperatures (T ~ 10 ® K) and low densities (jig ~ 10 “^ cm~^), typical values found 
in the outer regions of large X-rays clusters (Markevitch et ah 1996 ), the time-scale for 
electrons and ions to reach thermodynamical equilibrium through Coulomb collisions is 
about fei — 4 X 10 ® yr, comparable to the Hubble time. 

In order to study these low density regions where strong departure from 
thermodynamical equilibrium is expected, it has been pointed out by Kang et ah 
(1994) that an Eulerian code is well suited. In this paper, we therefore present a new 
Eulerian code and then use it to model non equilibrium processes during the formation of 
Zel’dovich pancakes. 


Pancakes appeared for a long time as fundamental tools in Cosmology. Zel’dovich 
(1970) was the hrst to point out that sheet-like structures could form through gravitational 
instability. These pancakes was hrst motivated by the neutrino-dominated scenario of 
structures formation, where they form naturally. Although this scenario seems today in 
difficulties with observations, walls and hlaments are still observed in both observational 
and theoretical studies based on other more popular scenarios like Cold Dark Matter or 
Mixed Dark Matter ( |Cen fc Ostriker 1992| , [Peebles 1993|) . Therefore, pancake geometry is 
not only an idealized case for testing numerical codes, but is also cosmologically relevant. 


Hydrodynamics of pancake collapse have been studied by several authors ([Bond et al.\ 
[1984[ [Shapiro fc Struck-Marcell 1985[ , [Anninos fc Norman 1994[) with the baryon component, 
obeying the hydrodynamics equations and coupled to collisionless dark matter particles. 
Different numerical schemes for modeling the hydrodynamical equations were used. Bond 
et al. (1984) and Anninos and Norman (1994) used an Eulerian scheme, while Shapiro and 
Struck-Marcell (1985) used a Lagrangian scheme. All these studies were dedicated to the 
calculation of the cooled mass fraction formed via pancake collapse. On the contrary to the 
previous studies, we focus in this paper on non-equilibrium phenomena which are enhanced 
in the large scale, low density part of pancakes. We introduce three temperatures to describe 
the thermodynamic evolution of electrons, ions and neutral particles. Chemical evolution of 
the primordial Hydrogen-Helium gas is solved without assuming ionization-recombination 
equilibrium. Moreover, we also model electronic conduction with a flux-limited diffusion 
scheme. 
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This paper is organized as follows. In section 2 we present the basics equations 
governing the system and we describe our hydrodynamical code with validating tests. In 
section 3 we present the results of simulations concerning the formation of pancakes. We 
show that, in the general case, a complex structure forms, with a thermal wave escaping 
out of the shocked region together with an equipartition wave where electrons and ions 
temperatures can differ signihcantly. We dually discuss in section 4 the role of various 
parameters, such as the baryons density parameter Qb, the wavelength of the initial 
perturbation L and the pancake collapse epoch Oc. 


2. PHYSICS AND NUMERICAL METHODS 
2.1. Basic equations 

The equations are written in comoving coordinates, through the transformation 
r = a{t)x where a{t) is the expansion factor. We take here a as the time variable. We 
assume an Einstein-de Sitter universe, with zero cosmological constant. The velocity of 
Dark Matter particles and fluid elements are respectively v = and u = Dark 
matter particles satisfy the equations of motion 


Du 

Da 


2 - D/2 


V 


3D, 


( 1 ) 


D is the background density parameter, and H the Hubble constant (both parameters are 
time-dependent quantities). For Lagrangian fluid elements, one has to add the pressure 
term 


Du 2-D/2 3D„ , 1 V^P 

TT =-“ ~TTH - 

Da a 2a^ pB 

and to consider the continuity equation 


( 2 ) 


1 Dp 
p Da 


- + v, 

a 


u 


( 3 ) 


The total pressure P = Pe + Pi + Pn D the sum of the electrons, ions and neutral particles 
partial pressures. The gravitational potential satisfies the Poisson equation 



^ = 5 


P 


( 4 ) 
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where p = + Pb the total mass density and p is the mean backgronnd total mass 

density. Thronghont this paper, densities refer to proper physical qnantities and not to 
comoving qnantities. We consider six chemical species of nnmber densities of each species 
null nHiii nHei, nHeii and nueiii- The partial pressnres are then related to the kinetic 
temperatnres as Pe = riekTe, Pi = {uhii + nneii + nneiii) kTi and Pn = (uhi + nnei) kTn. 

We dehne the specihc volnme V = 1/pb^ occnpied by a nnit mass of baryons. The 
specihc internal energy for each thermodynamical specy {a = e, i, n) follows then the 
eqnation of state for a mono-atomic gas (7 = 5/3) 

= -^PaV (5) 

7 - 1 

and satishes the hrst law of thermodynamics 


P>Ea 

Da 


\a J Da 


( 6 ) 


The hrst term in the right hand-side of eqnation (|^) is the PdV work dne to expansion 
and comoving compression and the last term is the net heat sonrce per nnit mass dne to 
different irreversible, non adiabatic processes. 


2.2. Thermodynamical processes 

The thermodynamical evolntion of the plasma is treated in a self-consistent way with 
its chemical evolntion, withont assnming ionization - recombination eqnilibrinm. The 
thermodynamical processes modeled in onr code are: shock heating for ions and nentral 
particles, cooling and thermal condnction for electrons and eqnipartition between electrons, 
ions and nentral particles. 

We have nsed the collisional ionization rate and the radiative recombination rate given 
by Cen (1992) which inclnde correction terms for very high temperatnres. As we have to 
deal with three different kinetic temperatnres for the gas, the actnal rates are obtained by 
nsing the reduced temperature of the two reactants (Draine 1980; Draine & Katz 1986). In 
the case of proton-electron interaction, this writes T^i = f nPIkltuPIiY practice, if Tg ~ Tj, 

\ 'IPlp~r'Ple / 

this implies a small correction of the order of me/rrip. Bnt in some extreme cases where 
Ti ^ Te, this correction is not negligible. We also considered classical cooling processes snch 
as ionization, recombination and line cooling, together with bremsstrahlnng and Compton 
cooling by the Cosmic Backgronnd Radiation. The cooling rates are again those of Cen 
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(1992), modified using the reduced temperature. Note that cooling results in an internal 
energy loss for the electrons only. 

Massive particles (ions and neutral) share most of the entropy deposition due to shock 
heating. The electrons are mainly heated by energy exchange with the latter species. 

As a matter of fact, in a perfect fluid, shocks are discontinuities in the flow, obeying 
Rankine-Hugoniot relations. They imply that the post-shock temperature for a given 
particle specy “i” is T oc rriiD’^, where T) is the up-wind fluid velocity in the rest-frame of 
the shock front. Consequently, the post-shock electron temperature is irte/mp, much lower 
than the ions post-shock temperature and negligible compared to the hnal equilibrium 
temperature. We therefore neglect electrons shock heating. Shock heating will be treated 
in the code using the artihcial viscosity method (Richtmyer & Von Neuman 1974), which 
includes a linear and a quadratic viscous term. This can be written in one dimension, and 
in one dimension only, as a viscous pressure, that we add to the usual ions (resp. neutral 
particles) thermal pressure. 


Pi,vise = Pi (Cie 6 * 26 ^) where e 


aAx 


(3 -|- aVa; • u) 


(7) 


The equivalent energy source term entering equation (13) for ions is given by 


DQ 


I,Vise 


Da 


Di,visc^ ( ^ 


( 8 ) 


Similar equations apply for neutral particles. C\ and C 2 are two constants determined a 
posteriori by numerical tests, and Ax is the mesh spacing. 

Electrons are heated by ions through Coulomb interactions, and by neutral particles 
through short-range forces. The equipartition rates are computed using the momentum 
transfer cross section of the different interacting species. For example, the net kinetic 
energy transfer rate per unit mass between electrons and protons is (Spitzer 1962) 


. ,rj, /^4(27r)^/2g4^V2jj^^^^N^ 

Da Da m, (fcTe.)^/^ ] PsaH 

where Tgj and InAep are respectively the reduced temperature and the Coulomb logarithm 
of the two interacting particles. The heat transfer rates between the other chemical species 
can be expressed in a similar form (Draine 1980; Chieze, Pineau des Forets & Flower 1998). 


When electronic temperature gradients are present in the flow, a net heat flux appears. 
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written here in its “classical” form 


Qcl = --KeV^Tg (10) 

a 

with conductivity coefficient (Spitzer 1962) 

ji5/2 

Ke = 1.84 X lO-*^ (11) 

In the case of very high fluxes (very steep gradients or very high temperatures and low 
densities), the flux saturates to a value corresponding to a free transport of the electron 
internal energy at a fraction of the electrons sound speed. We choose the flux-limited 
diffusion scheme described in Cowie and McKee (1977). The formulae we use here are the 
followings 


Qe 


Qd 

1 + qd/qsat 


where the maximum (saturated) value of the heat flux is given by 


( 12 ) 


/okT \ 

qsat = 0A{ -^ UekT, 

\7imeJ 


(13) 


The heat source (or sink) which enter equation (j^), due to electronic conduction, is Anally 
given by 


DQe 

Da 


1 


Vx • Qe 


(14) 


Ions are also able to transfer heat through ions conductivity, but the conduction coefficient 
is reduced by the factor (rue/mp)^^^, so this extra heat flux is much less effective. Moreover 
the average kinetic velocity for the ion gas after the shock front is lower than the shock 
velocity. Consequently ions do not cross the shock front, except for a few supra-thermal 
particles. The thermal flux due to ions is therefore neglected. In the opposite, the electrons 
sound speed is {mp/rrie)^^'^ times larger than the shock wave velocity 
|1966|) . The electron heat flux therefore crosses easily the compression front and pre-heats 
efficiently the cold gas ahead of the shock. Finally, ions will in turn be heated through 
Coulomb energy exchange with electrons. 


(Zebdovich fc Raizer 
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The capacity of electron to transport heat via conduction is dramatically limited in 
presence of magnetic helds. Since no evidence for strong magnetic helds has come to our 
knowledge, we expect electronic conduction to occur at this rate. But during the collapse, 
and more specihcally during shock wave formation, a strong enough magnetic held could be 
generated and reduce the effect of conduction. Therefore, we made simulations with and 
without electronic conduction. 


2.3. Numerical technics 


We present here our ID hydrodynamical Eulerian code. The extension of this code to 
a fully 3D hydrodynamical scheme is presented in a companion paper (IChieze, Alimi fc 
P?eyssier (1998)|) . 


Our code is based on the operator splitting method with four consecutive steps. The 
hrst step is called the Gravity step. It solves the Vlasov-Poisson equations for Dark 
Matter particles and calculates the gravitational potential. The second step is called the 
Lagrangian step, and solves the adiabatic Hydrodynamics equations in their Lagrangian 
form. The third step is called Eulerian step, it calculates the projected hydrodynamics 
quantities on the hxed Eulerian grid from the perturbed Lagrangian grid. The last step 
is called the Dissipative step; it computes all local dissipative processes, cell by cell, 
using the densities resulting from the two previous steps. Our code in its hnal version is of 
second order accuracy both in time and in space. It allows great stability and efficiency. 
We present here its general features. We then show tests which demonstrate its ability to 
handle cosmological simulations. 


2.4. General Presentation of the Code 

We consider a two-fluid system (Dark-Matter and Baryons). The physical variables 
associated to Dark Matter particles (superscript j) are , the position of particle j, , 
the velocity of particle j. The discrete values of the flow on the grid (superscript i) are 
M\ the total baryons mass in cell i, Nx, the total numbers of particle X in cell i, SI, S'*, 
S'*, respectively the total entropy of electrons, ions and neutral particles in cell i, wj,, the 
velocity in the x-direction of interface i and hnally rj,, the position in the x-direction of 
interface i. 


Mass, particles numbers and entropies are zone-centered, while the velocity is 
face-centered. This is the well-known staggered mesh method ([Stone fc Norman 1992|). It 
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allows better accuracy when computing finite differences, and also reduces the number of 
interpolation to calculate fluxes, which are defined at cell interfaces. The variable is 
usual in pure Lagrangian schemes, it allows to compute densities, while the mass remains 
constant. As we will see below, it also improves the accuracy of both time integration and 
flux interpolation. 

The specific entropy Sa is defined for each thermodynamical specy a as Sa 
Using the energy equation (^, the time-derivative of Sa reduces then to 

_ SgBQa 
Da Eg Da 

The thermodynamical evolution of baryons is computed using the entropy 
(^) rather than the energy equation (^. This method does not introduce any 
spurious dissipation effects due to expansion or compression. Consequently, an 
flow ( BQg/Ba = 0 ) is strictly adiabatic during the Lagrangian step. 


= EgV^-^. 


(15) 

equation 

numerical 

adiabatic 


2.5. Gravity Step 


We evolve the Dark Matter particles using a Particle-Mesh algorithm ([Hockney fc 


Eastwood 198 Ij) with a Predictor-Corrector scheme. This time-integrator ensures both 
second-order time-accuracy and variable time-stepping. Let us suppose that particles 
positions and velocities are known at a given time a. We compute the predicted positions 
at time a + Aa with + Aa uL The superscript (1) means that the quantity 

is the predicted quantity evaluated at time a + Aa. The Dark Matter density fields p* 
and are then evaluated with a CIC interpolation scheme. In order to solve now 

the Poisson equation, we need to know also the predicted Baryon density field (p^)^^^- 
This one is deduced from the Baryon density and velocity fields at time a by solving the 
continuity equation. From the gravitational potential (j) at time a and from its predicted 
quantity at time a + Aa, we deduce then the forces E^ and for the Dark 

Matter particles by inverse interpolation, at respectively positions x^ and positions 
Velocities and positions for Dark Matter particles are finally updated according to the 
formula 


(.yi)(2) _ yj 


2 - D/2 (F)(2) +t;i 


3D E^ + (F^')(i) 


Aa 


a 


2 


2 


( 16 ) 
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Aa 2 


(17) 


where a is now the centered expansion factor a + Aa/2, is the corresponding density 
parameter and and are the updated position and velocity for particle j at 

time a + Aa. 


2.6. Lagrangian Step 

At this step we solve the adiabatic Hydrodynamics equations, together with shock 
heating and electronic conduction. All the quantities involved in these equations depend 
on spatial derivatives of the flow. These ones are estimated with a hnite-difference scheme. 
Chemical, cooling and equipartition processes which are purely local (i.e. they don’t depend 
on spatial derivatives of the flow) are solved in the Dissipative step. 

The adiabatic Hydrodynamics equations for interface i write for an explicit scheme 
(hrst order: superscript (1)) as 


(^) 



(18) 


_ 2-D/2^. 2Ax^P^-P^-^ 

“ a A^ aH^ M* + M*-i 

where the pressure P* includes electron, ions and neutral partial pressures, and the artihcial 
viscous pressures for shock heating (equation |^). The entropy equations due to electronic 
conduction for electrons, and to shocks heating for ions and neutral particles for the cell i 
are derived in a similar way by hnite differencing equations (^, (^) and (|T^. Only the 
variables (r^, u^., SI, SI, S\, i = 1,N) are coupled, the total mass and total numbers of 
particles remains strictly constant during the Lagrangian step. In order to integrate the 
previous equations, we have to compute the increments (designed by A) of each variables 
between time a and time a + Aa. Several methods are possible. The most straightforward 
one is to use directly the explicit estimation we already mentioned (hrst order method; see 
Stone &: Morman 1992| ). It however needs a rather strong condition on the time step in 
order to be stable, namely the Courant condition. The implicit method, on the other hand, 
allows much larger time-steps. It is very stable, but very CPU-time consuming, because 
it implies to invert a band matrix. In this paper, we have preferred to use a second-order 
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time integrator scheme inspired by the implicit method. It consists to derive second order 
increments by Taylor expanding the hrst order increments given by equations (|^ and 
(0)> with respect to the flow variables, namely Se, St and S'„. For example , the 

acceleration terms for interface z, which depends on u^., r*, S'*, 

S'*, S'^, and yields the following second order (superscript (2)) estimation of the 
velocity increment 


Aa 


( 2 ) 


f Am, 


( 1 ) 


\ Aa 


+ 


Aa 

~Y 

Aa 


' Ar, 


i-i\ (1) 


d ( Am,.*' 


( 1 )' 


\ Aa ) dr^~^ y Aa J 

/Am,*-^\^^^ d fAujV^^' 


y Aa 


du^~^ y Aa 


( 20 ) 


where only the partial derivative with respect to the and m^“^ terms has been written 
for sake of simplicity. It is important to add partial derivatives of the first order increment 
with respect to all variables involved in equation (P^). This method differs from the 
standard “predictor-corrector” scheme, since the second order correction presented here is 
computed analytically using partial derivatives of the first order increments. At the end of 
the Lagrangian step, we have the new entropies, the new velocities and the new interface 
positions. Mass and numbers of particles have not been modified. 


Finally the time-step is controlled using usual methods (see e.g. [Stone fc Norman 


1992[), with a constraint for each of the following processes: artificial viscosity, conduction. 


gravity and gas dynamics. For example, for a pure gas dynamics problem, the time-step is 
controlled by the following criterion 


Aa 




^^HAx 


0 mm 


'{a^HuY -I- -I- 


( 21 ) 


where Cg is the adiabatic sound speed, and c^isc is related to the viscous pressure by the 
formula 


^visc 



( 22 ) 


The factor 4 in the last equation corresponds to the stability criterion for a viscous fluid 
with constant viscosity coefficient u] the equations of motion in this case are analogous to 
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the well documented diffusion equation, for which the stability criterion is established as 
At < (Aa;)^/(4z/) (see Stone & Norman 1992). The Courant safety coefficient Cq has to 
be chosen less than 0.5 in order for an explicit scheme (first order) to be stable, and even 
smaller (typically 0.1) to be accurate. On the contrary, our second order scheme remains 
stable and accurate, even with the ultimately large time-stepping given by Co = 1 (see the 
tests section). 


2.7. Eulerian Step 

We now want to re-map all variables from the disturbed Lagrangian grid to the hxed 
Eulerian one. This step has to conserve mass, numbers of particles, internal energies 
of electrons, ions and neutral particles, and momentum. We dehned the left-centered 
momentum in cell i as PI = M^u], and the right-centered momentum in cell i as 

in order to deal only with zone-centered quantities. The projection step 
consists then to solve the advection equation (written here only for the mass) 

^ J pdV = — j pu ■ dS (23) 

for each zone which has a control volume E* = (Ax)^ ~ '^x)- A hnite difference 

approximation of this integral equation is 

- M* = AM* - AM*+^ (24) 

where Mp are the projected masses on the Eulerian grid and AM* is the advected mass 
through interface i. This scheme is strictly conservative for the projected variables by 
construction. The main problem of the advection procedure arises when one calculates the 
total mass contained in the advected volumes. To do so, we have to calculate a realistic mass 
distribution within each cell, and then integrate this distribution in the advected volume. 
There are basically three tractable methods: uniform, linear and parabolic distributions. 
These three methods are known as the Donor Cell method, the Van Leer method (Van 
Leer 1977) and the Piecewise Parabolic Interpolation (PPI) method (Woodward & Colella 
1985). The first one is very simple but quite diffusive. It is first order accurate in space. 
The two other methods are much less diffusive and are respectively second and third order 
accurate in space. The PPI method is the most time consuming, while Van Leer method 
offers a good compromise between accuracy and efficiency (see Stone & Norman (1992) for 
a general description). 
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We calculate interpolation functions for mass, momentum and total internal energy 
{E = Ef. + Ei + En) only. The numbers of particles within the cell are distributed according 
to the mass distribution, and internal energies for the 3 species are distributed according 
to the total internal energy distribution. This ensures exact mass and charge conservation 
within each cell, and avoids spurious decoupling between the 3 temperatures. The entropy 
are then updated from the projected internal energy, and the velocity is updated from 
the left-centered projected momentum, the right-centered projected momentum and the 
projected mass according to the formula 


pi—1 


+ Pl 


-F M* 


(25) 


This equation yields an exact momentum conservation. All new Hydrodynamics variables 
on the Eulerian grid are now known. 


2.8. Dissipative Step 


We solve now cell by cell, all purely local dissipative processes: chemical reactions, 
thermo-chemical energy exchanges, equipartition and cooling. All these processes are 
described by very stiff equations. They imply necessarily very short time-steps, which would 
slow down dramatically the whole simulation and increase the CPU-time up to prohibiting 
values. Consequently to avoid this time-step catastrophe, we solve these stiff equations by 
using n consecutive sub-steps. At each sub-step and for each cell, we invert a 6 x 6 matrix 
for chemical reactions of the 6 involved species and a 3 x 3 matrix for the entropies. Our 
algorithm is fully vectorized and very effective, allowing high accuracy and stability. The 
number of sub-steps depends on the physical state of the cell. A similar method was used 
by Anninos and Norman (1994). Such a method is not justihed as soon as the time scale of 
dissipative processes is much shorter than the dynamical time scale given by the Lagrangian 
step. This happens in very dense regions where cooling can be then overestimated. We 
discuss how to avoid such a mistake in [Chieze, Alimi fc Teyssier (1998) . 


In the Eulerian 


low resolution case presented here, this happens only in one or two cells at the pancake 
center. However, using the prescription presented in [Chieze, Alimi &: Teyssier (1998)| (which 
consists essentially to turn off line cooling in the central cell only), we hnd that cooling at 
small scales has little influence on the large scale flow. 


2.9. Numerical Tests 
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2.9.1. Advection Test 

This test was proposed by Stone & Norman (1992) to qualify the advection scheme 
only. We consider a box of length unity hlled with a gas of homogeneous density. With 
a resolution of N = 512 cells and periodic boundary conditions, we model the advection 
at constant velocity of a single square pulse of density 10, initially sampled by 50 cells. 
Our results obtained with the three different projection schemes are very close to those 
obtained by Stone & Norman (1992) at the time when the pulse has crossed half of the 
computational space. The Donor cell interpolation is dramatically diffusive, while the two 
other schemes reproduce the sharp features relatively well. 


2.9.2. Shock Tube Test 

This test, also called the classical Riemann problem, is may be the most widely used to 
qualify hydrodynamics codes (Sod 1978). The initial conditions we use here are similar to 
those used by Stone & Norman (1992). We consider a box of size L = 1, N = 128 cells and 
u = 0. We separate the box in two regions (left and right) with the following conditions 
Pl = 1, Pl = I and pr = 0.125, pr = 0.1. We assume 7 = 1.4 and reflective boundary 
conditions. We use for this test the quadratic term of the artihcial viscosity only, with 
Cl = 0 and C 2 = 3. In hgure we plot the different prohles obtained at time t = 0.245 
for our standard simulation parameters, namely a Courant safety factor Cq = 0.5, the Van 
Leer advection scheme and our second order time integrator. The results of our code is 
comparable to other methods, and very close to the analytical solution. Note however 
that the specihc energy in the post-shock region is slightly better reproduced here than for 
example in Stone & Norman (1992). We think that this better agreement is mainly due to 
our second order time integrator. 

To analyze the effect of our second order time integrator, we run the same simulation, 
but with a Courant safety factor Co = 1, which is the ultimate possibility. We plot in 
hgure (^ the specihc energy prohles obtained at time t = 0.245 for the three diherent 
advection schemes (Donor Cell, Van Leer and Piecewise Parabolic) using the second order 
time integrator, together with the prohle obtained with the Van Leer advection scheme, but 
using the hrst order, explicit method. Note that in this latter case the solution is strongly 
unstable, while for the three former cases, the solution is identical to the Co = 0.5 case. 
This illustrates the interest of our second order time integrator. We also learn from this 
hgure that the Donor Cell advection scheme is dramatically dihusive, and is practically of 
no use. The improvement of the solution from the Van Leer to the Piecewise Parabolic 
scheme is real, but not very dramatic, although for the latter the computational cost is 
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much higher. This justify the use of the Van Leer method as a standard choice. 


2.9.3. Blast Waves Test 

This test was used by Woodward & Colella (1985) to compare different hydrodynamics 
codes. It is may be the most adapted test for cosmological applications, because strong 
shocks are generated with Mach numbers ~ 10^ and are interacting violently. The initial 
conditions we use are the following (see Stone & Norman 1992) L = 1, N = 1200 cells, p = 1 
and u = 0. In the left tenth of the box, we put p = 1000, in the right tenth p = 100 and 
in the middle p = 0.01. We assume 7 = 1.4 and reflective boundary conditions. To have a 
more precise description of the different features involved in this simulation, see Woodward 
& Colella (1985). For this test, the linear term of the artihcial viscosity is necessary to damp 
small oscillations occurring otherwise after strong shock fronts (we assume Ci = 1 and 
C 3 = 3. We consider a Courant safety factor Cq = 0.5 and our second order time integration 
scheme. Had we decided to use here the hrst order time integrator (explicit method), the 
Courant safety factor would have been Co — 0.1 in order to recover similar results. We 
show in hgure (j^) the density prohles obtained at time t = 0.038 in order to compare with 
Stone & Norman (1992) and Woodward & Colella (1985). At that time, the two shock 
waves have already interacted at x ~ 0.7 and are moving back to their original positions. 
The Piecewise Parabolic scheme presents the best agreement with Woodward Sz Colella 
(1985) reference simulation, although Van Leer results appears to be very similar too. As 
already mentioned. Donor cell scheme is unable to handle correctly sharp discontinuities, 
but more dramatic is the bad positioning of the shock fronts. During this run, total energy 
conservation was for Donor cell. Van Leer and PPI schemes respectively 7.6%, 2.1% and 
1 . 1 %. 


3. PANCAKES COLLAPSE 
3.1. Initial Conditions 

We take for all runs the same initial conditions, beginning at epoch Zi = 200 where 
matter and radiation are well decoupled. The total density contrast has the following 
distribution 


^(a:) = ^cos (27ry) 
a. \ LJ 


( 26 ) 
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and the baryons density distribntion is taken eqnal to Pb{x) = ps {I + S{x)). Dark matter 
particles are then moved according to Zel’dovich approximation, with the displacement held 
corresponding to the above density contrast. Initial velocity held for baryons is given by 
the linear growing mode 


Ux{x) = -(27) 

a 

The temperatnre for ions, electrons and nentral particles are chosen nniform and 
eqnal to = Tq (1 + Zi), where Tq = 2.7 K is the cosmic radiation backgronnd 

temperatnre today. We snppose that the helinm mass fraction is D = 0.24. The initial 
ionization fraction are taken from Peebles (1993). In all onr rnns, we take D = 1 and 
Ho = 50 km Mpc~^. Therefore, the single cosmological parameter of interest here 
is Qb- We have adopted periodic bonndary conditions to ensnre total mass conservation. 
We nsed Van Leer projection scheme since it is a good compromise between accuracy and 
CPU time and a Courant safety factor Cq = 0.5. We consider in the following a reference 
pancake dehned by a comoving length L = 16 Mpc h~^, a collapse epoch Oc = 0.2 and a 
baryon density parameter Qb = 0.1. 

3.2. Adiabatic Collapse 

In this standard case, the gas is assumed to be fully ionized, and described by a single 
kinetic temperature. Sunayev & Zel’dovich (1972) and Bond at ah (1984) have derived 
analytically density and temperature prohles in that case. They have shown that a shock 
wave appears off center at a very small radius (typically a few 10“® Mpc) which corresponds 
to the sonic radius. The flow is almost hydrostatic in this inner, unshocked region of high 
density contrast. The accretion shock front propagates outward, leaving an almost uniform 
pressure all over the accreted gas. Here, contrary to Bond et al. (1984) and Shapiro & 
Struck-Marcell (1985), we don’t resolve this very central region, because we are using a 
regular mesh (this mesh has been chosen in such a way to describe better the large scale, 
outer region of the pancake). We are aware of the fact that the better is the resolution, 
the higher is the density contrast in the central cell, until the resolution corresponding 
to the sonic radius is reached. For adiabatic runs, this has no incidence on the general 
aspect of the flow, since the pressure is almost uniform, in contrast with non adiabatic 
runs. We present in hgure (|) velocity, temperature, pressure and baryons density prohles 
at z = 0 for various mesh resolution. It can be seen, besides an increasing sharpness of 
the shock front and an increasing density contrast in the central cell, that the resolution 
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has little quantitative influence on the results. The velocity held is typical of an accreting 
quasi-hydrostatic flow. The usual self-similar scaling laws {Tad cx: Uad oc and 

Pad — const) are well reproduced. 


3.3. Non Adiabatic Collapse 

We turn now to the analysis of non adiabatic pancake collapse, focusing on the 
temperature structure in the flow resulting from the microscopic collisional processes among 
the various species. We examine both the effects of energy exchange processes and of 
electronic conduction. 


3.3.1. No Electronic Conduction 

We follow precisely chemical reactions with the corresponding thermochemistry, 
collisional cooling processes and equipartition processes, as discussed in Section 2. We 
first suppose that electronic conduction is ineffective. Figure @ shows that the electronic 
temperature decouples from the ion temperature at 600 kpc h~^ from the mid-plane. The 
maximum departure from equilibrium is found near the shock front at 1.1 Mpc h~^. In this 
region, ion temperature is about 10^ K when electron temperature barely reaches 10® K. 
Furthermore, the ion and electron temperature profiles are different, with opposite gradients: 
the electron temperature steadily drops towards the front, while the ion temperature rises. 

Non equilibrium chemistry is required especially in the post-shocked regions in which 
ionization gradually reaches its near-equilibrium value. Recall that electrons are essentially 
heated behind the shock front by Coulomb collisions with hot ions. Since equipartition 
processes are conservative, the total energy density (i.e. the total pressure) is unchanged 
relative to the single temperature case. By themselves, equipartition processes do not affect 
the dynamics of the flow, which is only modified - relative to the adiabatic case - by cooling. 
Since the equipartition rates per particle are proportional to density, the temperature of 
ions and electrons are well coupled in the dense parts of the pancake, but significantly 
decouples in the low density outer regions. The equilibrium point (where Tg and Tj differ 
from less than 2%) will be roughly estimated by analytical calculations presented below. 
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3.3.2. Electronic Conduction 

In that case, we suppose that electronic conduction is fully effective, with flux-liniited 
diffusion (Cowie & Me Kee 1977). We hnd that conduction is saturated only in the central 
region, where temperature gradients are very stiff, and in the most outer regions, where the 
density is very low. Computing thermal fluxes is a complicated task, since it depends on the 
ionization state of the gas, which in turn depends on the propagation speed of the ionizing 
thermal wave. We need therefore a good sampling of the ionization front to accurately 
track the thermal front. This explains why we decided to choose a linear mesh, in order to 
achieve a fair sampling of temperature and abundance prohles in the thermal wave. 

We plot in hgure the pancake state at z = 0. Note that the ionization state of the 
gas is very well sampled and that abundances gradually evolve towards their equilibrium 
values. The shock front position at 1 Mpe h~^ is very close to its position in the “no 
conduction” run. The main effect of conduction is the thermal precursor, which pre-heats 
and pre-ionizes the gas up to 1.5 Mpe h~^ from the mid-plane. Because equipartition 
processes are slow in this region, ion temperature reaches only 10^ K, while electron 
temperature is 10® K. Nevertheless, this results in a slight dynamical effect on the flow: 
ion pressure gradients cause a small deviation from the pressureless velocity prohle in the 
unshocked region, clearly visible on this hgure. Shock-heated regions occupy a total volume 
of 2 Mpe h~^, while unshocked but pre-heated (Tj ~ 10® K) and pre-ionized (Te ~ 10® K) 
regions occupy a total volume of 1 Mpe h~^. As we will see in the next part, the efficiency 
of electronic conduction depends strongly on the pancake size. 


3.4. Varying the Pancake Parameters 

In this section, we study the inhuence of the different parameters on the pancake 
structure. We develop approximate formulae which guide us for our conclusions. We assume 
that, in any case, the dynamical state of the pancake is given by the adiabatic collapse. 
Following Bond et al. (1984) and Shapiro & Struck-Marcell (1985), it is then possible to 
derive interesting formulae for the pancake evolution. 

We assume hrst that each huid element q is shock-heated at an epoch given by 

^ _ 7rg 
ttc sin Tiq 

This corresponds to the epoch when the corresponding collisionless particle crosses the 
center of the pancake. We suppose also that the gas is fully ionized and that the flow 


(28) 
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remains strictly adiabatic. Before the shock front, the flow follows the pressureless solution 
of pancake collapse (|Zehdovich 19701 ) . Using Rankine-Hugoniot relations and assuming 
that the post-shock peculiar velocity vanishes, we get the post-shock temperature and the 
post-shock over density 


kT = —larupa 

(29) 

1 ^ ^ cos7rg\ 

I -I- (5 = 4/ I — vrg- 

\ sin 7rg ) 

(30) 


As we assumed that the peculiar velocity is zero in the post-shock region, the temperature 
evolves afterwards as a~^ and the density as a~^. Mass conservation implies also that the 
Eulerian position of a given fluid element in the shocked region is given by 


/’'? / cosvrg 

x= 1 — Trg- 

JO \ sm vrg 


(31) 


This allows us to describe the dynamical evolution (single temperature case) of the 
pancake. Let us now estimate the thermodynamical evolution using our three temperatures 
formalism. The equipartition time-scale for electron-ion energy exchange is given by 


rj^3f‘2 

4 - = 503 f , sec. (32) 

He In A 

Assuming that the effective temperature Tg* is equal to the post-shock adiabatic 
temperature, the equipartition time-scale remains then constant during the post-shock 
evolution (we neglect the slow variation of the Coulomb logarithm, taking InA ~ 40). It is 
now possible to solve analytically the equipartition equation 


4(Te-T,)+2-(Tg-T,) = -A(Tg-T,) 

at a tg^i 

The single temperature, given by equation (^Of) , is related to Tg and Tj by 

^ _ UgTg -F njTj 
Tig -t- Hi 


(33) 


(34) 


We hnally obtain semi-analytical temperatures prohles, plotted in hgure (^), for the 
reference pancake at z = 0. Note that numerical results agree qualitatively well with 
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our analytical theory. Because we assumed that the post-shock velocity vanishes, we 
overestimate the shock front position and the post-shock temperature. In the numerical 
calculation, Lagrangian fluid elements pile up deeper than we assume in our analytical 
calculation. This explain the visible differences between numerical and analytical results. 

The pancake structure is fully described by three characteristic points, as we already 
presented in the upper sections: the compression point where ions and neutral particles are 
shocked, the thermal point which marks the end of the thermal wave and the equilibrium 
point where equilibrium between electrons and ions is recovered. In the small q limit. Bond 
et ah (1984) derived the Lagrangian coordinate of the compression point 


qs oc 



(36) 


We can write, assuming that a — Oc ^ ac 


qs oc 


a — Qc 

dc 


1/2 


(36) 


We derive here (for g <C 1 and a — Oc ^ Oc) the Lagrangian coordinate of the equilibrium 
point (where electrons and ions temperatures differ from only 2%) 



It appears that this coordinate depends strongly on the pancake size. This is due to 
the strong dependence on temperature (T^/^) of the equipartition time-scale. When the 
baryon density parameter is decreased from 0.1 to 0.01, the equilibrium point reaches much 
deeper regions. The hnal mass at thermodynamical equilibrium is about 85% of the total 
shock-heated mass for the VLb = 0.1 case, and decreases to about 50% for the VLb = 0.01 
case. The time-dependence of the equilibrium point is much slower than for the shock front. 
This means that non-equilibrium features appear mainly at late epoch and in the most 
outer regions where density has decreased (and temperature has increased) sufficiently. 

Formula (^) is valid only for the small q regime. To test the validity of our description 
at later epoch (the large q regime), we run numerical calculations with different pancake 
sizes. All these runs were computed assuming that electronic conduction is effective, and 
that VLb = 0.1, Oc = 0.2. We plot in figure (j^) the shock front Eulerian position Xg, 
the equilibrium point position Xgg and the thermal front position Xth obtained in these 
simulations. The shock front position is almost independent of the pancake size, as expected 
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by formula (^). The equilibrium point is found in deeper regions as L increases, as 
expected by formula (|37[). The thermal wave has approximately the symmetrical behavior 
than the equipartition wave. 

For L ~ 10 Mpc h~^, cooling starts to be important and the flow does not remain 
adiabatic. The main consequence is that the compression point propagates outwardly at a 
lower velocity than for the adiabatic case. This characteristic pancake length corresponds 
roughly to pancakes for which the average cooling time is equal to the Hubble time. Indeed, 
assuming that the mean temperature is given by equation with g = 1/4, 5 ~ 10 and 
using the equilibrium cooling curve to calculate the cooling rate, we hnd Lcooi — 8 Mpc h~^. 


4. CONCLUSIONS 

In this paper, we use a three temperature formalism to describe the thermodynamical 
evolution of pancakes. We show that the assumption of thermodynamical equilibrium 
between ions and electrons is not valid in the outer regions of pancakes. The corresponding 
temperature prohles show differences up to one order of magnitude near the shock front 
and reach thermodynamical equilibrium near the centre, where the overdensity is 5 ~ 10. 
Electrons and ions decoupling is stronger for large pancake sizes or for low VLb- This could 
have observational consequences for X-ray temperature prohles. The dynamical (total) 
pressure in the outer regions of clusters can differ from the observed electrons pressure up 
to one order of magnitude. This effect can be even stronger for non-relaxed clusters, where 
recent mergers strongly decoupled electrons and ions through shocks. This could lead to an 
underestimate of the cluster mass in this regions. In the central hydrostatic part of clusters, 
where thermodynamical equilibrium is efficiently recovered, such effect are not likely to 
appear. We address this question using a fully 3D X-ray clusters modeling in a companion 
paper (Chieze, Alimi & Teyssier 1998) 

In the outermost regions, we show that electronic conduction, if effective, leads to a 
thermal wave escaping the shock-heated regions of the pancake. This thermal precursor 
could have interesting cosmological consequences, such as late reionization and heating of 
non collapsed regions. In our case, the precursor is strongly conhned by the very high infall 
velocity towards the pancake centre. In a more complex geometry, one can imagine that 
the size of the precursor could be much more extended, leading to efficient heating of the 
intergalactic medium. This effect might be detected in the vicinity of very large X-rays 
clusters. 


Non-equilibrium regions (from the thermal precursor to the equilibrium point) are 
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very extended and dominate in volume the structure of pancakes (this paper) and X-ray 
clusters (Chieze, Alimi & Teyssier 1998). Using high sensitivity experiments, like the X-ray 
satellite XMM (Fabian 1987), we will be able in the future to observe the low density part 
of clusters, where all these processes are likely to be very important. 

The authors would like to thank an anonymous referee for the constructive remarks 
that allow us to increase the quality of our work. 
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Fig. 1.— Density, pressure, velocity and specific energy profiles obtained for the shock tube 
test at time t = 0.245. We use for that run our standard parameters: Van Leer advection 
scheme, second order time integrator, Courant safety factor Cq = 0.5. The analytic prohles 
are shown as solid lines. 
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Fig. 2.— Specific energy profiles obtained with the shock tube test for different advection 
schemes, second order time integration and the “ultimate” Courant safety factor Cq = 1. We 
also plot the solution obtained with the explicit method (first order) for the same Courant 
safety factor. The analytic prohle is shown as the solid line in each graph. 












collapse of the reference pancake (L = 16 Mpc Oc = 0.2 and his = 0.1. In each graph 
we plot three runs with increasing resolution {N =128: dashed line, 256; dotted line and 
1024: solid line). 
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temperatures, ionization fraction for HI (dashed line), Hel (dotted line), and Hell (solid 
line) at z = 0 for the reference run without electronic conduction. 
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Same as figure (5) for the reference run with electronic conduction. 
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Fig. 7.— Analytical ions (dashed line) and electrons (solid line) temperatures profiles for 
the reference pancake at z = 0. Compare this hgure to the numerical results obtained in 
hgure (5). 








describe non-equilibrium features in the flow for various pancake sizes. 





